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1.  INTRODUCTION 


Modem  scientific  and  engineering  projects  have  become  too  complex  and  costly  to  be  carried  out  with 
just  the  old-fashioned  "trial  and  error"  approach.  Therefore,  it  is  imperative  that  such  projects  first  be 
modeled  and  verified  mathematically.  Unfortunately,  due  to  the  complexity  of  these  projects,  the 
mathematical  modeling  seldom  can  be  done  entirely  in  analytical  forms.  As  a  consequence,  one  must 
resort  to  "approximate"  methods;  that  is,  methods  which  either  entirely  or  partially  use  numerical  methods 
when  modeling  the  project.  The  most  popular  approximate  methods  in  use  today  are  the  finite  element 
and  the  finite  difference  methods.  A  new  method  is  the  finite-volume  method,  however,  it  is  not  yet  very 
popular  in  modeling;  this  is  very  likely  due  to  the  fact  that,  as  the  models  go,  it  is  rather  laborious  to  use. 
Consequently,  in  this  brief  description,  attention  will  be  paid  more  to  finite  element  and  finite  difference 
methods  than  to  the  finite  volume  method.  The  emphasis  here  is  to  give  an  overview  of  these  methods 
so  that  interested  readers  can  pursue  details  in  the  literature  on  their  own. 

In  section  2,  some  basic  notions  and  applications  of  finite  element  and  finite  difference  methods  are 
given.  Here  also  the  finite  volume  method  is  briefly  described.  Merits  of  one  vs.  the  other  method  are 
discussed  in  section  3.  Section  4  is  mainly  devoted  to  a  listing  of  possible  potential  applications  of  these 
methods,  and  contains  the  discussion  and  conclusion.  In  the  appendix,  listings  of  the  available  literature 
on  these  methods  are  given. 

2.  EXPOSITION  OF  THE  METHODS 

Here  simple  applications  of  finite  element  and  finite  difference  methods  are  described,  while  a  brief 
overview  of  the  finite  volume  method  is  given. 

2. 1  Finite  Element  Method.  A  mathematical  model  of  a  physical  system  normally  involves  a  number 
of  variables  and  functions  fex(x)  representing  fields,  velocities,  etc.  Here  %  represents  the  coordinates 
of  the  domain.  The  problem  is  that  the  exact  function  fgxOi)  is  not  known  for  all  the  points  x.  Hence, 
one  introduces  an  approximation  of  fex(^)>  denoted  simply  as  f(X).  The  error  function  eOl), 


e(x)  =  f(x)  -  fex(x). 


(2.1.1) 


1 


measures  the  quality  of  the  approximation.  To  construct  an  approximate  f  (x)  it  suffices  to:  (1)  write  an 
expression  containing  n  parameters  aj. 


f(x)  =  f(x;  aj,  ....  a^),  (2.1.2) 

and  (2)  relate  (determine)  these  parameters  to  n  values  of  fex(x)  in  the  domain,  which  may  be  known  or 
may  stiU  have  to  be  determined  by  some  other  methods  (Dhatt  and  Touzot  1984).  Formally,  this  may  be 
achieved  by  forcing  the  error  function  e  (X)  to  be  zero  at  "n"  points  in  the  domain. 

The  question  that  immediately  arises  is  how  to  construct  the  approximate  function  f  (%;  a^,  ...,  a„), 
dependent  on  parameters  aj,  i  =  1,  ...,  n.  Rather  frequently,  an  approximation  function  is  chosen  to  be 
a  linear  function  of  parameters  aj,  i  =  1,  ...,  n: 


^2 


fOt)  =  <Pi(lt)P2(x)-P„(j^)> 


1“ 


n 


J 


(2.1.3a) 


(2.1.3b) 


or  in  shorthand. 


f(X)  =  <  P(X)  >  (a„)  , 


(2.1.3c) 


where  PjOi),  i  =  1,  ....  n,  are  linearly  independent  complete  sets  of  functions.  In  the  finite  element 
method.  P’s  have  been  chosen  as  polynomials,  although  other  sets  of  functions  may  be  used.  Parameters 
aj,  i  =  1,  ...,  n,  are  the  generalized  parameters  of  approximation. 
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Unfortunately,  parameters  a^  generally  do  not  have  a  direct  physical  meaning.  Thus,  as  far  as  the  finite 
element  method  goes,  they  are  conveniently  replaced  with  nodal  values  of  the  function  fex(x)  at, 
points  with  coordinates  Xj,  —  •  ^m-  nodal  approximation  is  further  required  to  satisfy  the 

following  relations: 


f(Xi)  =  fg^CXj)  =  fi,  i  =  1,  2,  ...,  m. 
Hence,  the  approximate  function  f(x)  can  now  be  written  as 


f(x) 


Ni(x)fi 


ft 

h 


=  <Ni(x)N2(X)  ...  N^(x)> 


(2.1.4) 


(2.1.5a) 


(2.1iib) 


or,  in  shorthand  notation. 


f(x)  =<N(x)>(f^).  (2.1.5c) 

Here  fj,  i  =  1,  2,  ...,  m,  are  nodal  parameters,  and  functions  N(x)  are  called  (nodal)  interpolation 
functions.  Qearly,  consistent  with  relations  (2.1.4)  and  (2.1.5),  one  has 

Nj(Xi)  =5ij,  e(Xi)  =0.  (2.1.6,7) 
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Understandably,  for  some  approximate  function  the  generalized  approximation  (2.1.3)  and  the  nodal 
approximation  (2.1.5)  must  be  equal.  Hence  <  P(x)  >  and  <  N(x)  >  must  be  related  to  each  other. 
First  of  all  from  relation  (2.1.3)  one  can  write. 


Fi  =  <  P(x)  >  (a),  i  =  1,  2,  ...,  m. 
In  matrix  form  the  same  thing  is  written  as 


<  Pj  (Xj)  P2  (Xi)  ... 

PnO!l)> 

f  \ 

f2 

* 

<  Pi  (X2)  P2(X2)  - 

Pn(*2)  > 

h 

h 

• 

<Pl(xJP2(Xm)  - 

P„(Xm)  > 

) 

k  V 

(2.1,8a) 


(2.1.8b) 


where  one  should  notice  that  the  matrix  itself  is  not  generally  a  "square"  matrix.  In  a  shorthand  notation, 
this  relation  is  rewritten  as 


(f^)  =  [P]  (a„).  (2.1.8c) 

However,  unless  something  explicit  forbids  it,  one  may  choose  the  number  of  P-functions  and  the  number 
of  nodal  points  to  be  the  same. 

Relations  (2.1.8)  relate  nodal  parameters  f|,  i  =  1,  2,  ...,  m,  and  the  generalized  parameters 
aj,  j  =  1,  2,  ...,  n.  Substitution  of  (2.1.8c)  into  (2.1.5b),  after  comparison  with  (2.1.3),  yields 

<  P®  >  =  <  N(X)  >  [P],  (2.1.9) 
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and 


<  NOt)  >  =  <  P(X)  >  [Q],  [Q]  =  {PTK  (2.1.10a,b) 

Since  the  P-fiinctions  can  be  chosen  as  simple  polynomial  functions,  it  is  relation  (2.1.10a)  that  is  of  real 
value;  it  defines  the  interpolation  (nodal)  functions  in  terms  of  P’s.  The  problem,  however,  is  evaluating 
the  matrix  [Q]. 

The  following  is  an  example  of  the  construction  of  interpolation  functions  for  the  simplest  of  the 
elements,  the  linear  (two  nodes)  element;  it  is  linear  because  we  take  the  number  of  polynomial  basis 
functions  to  be  equal  to  the  number  of  nodes:  n  =  m  =  2.  The  one-dimensional  (1-D)  two-node  elements 
are  exhibited  in  Figure  1. 


X1 


X 


Figure  1.  The  1-D  two-node  element. 
The  polynomial  basis  functions  are  given  by  a  two-component  vector: 

<  P(x)  >  =  <  Pj (x)  P2 (x)  >  =  <  1 X  >. 


According  to  (2.1.8b  and  c)  the  matrix  [P]  is 


(2.1.11) 
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[P]  = 


1  Xi 

1  Xo 


This  matrix  can  be  easily  inverted  with  the  result 


i 


X2  -Xi 
-1  1 


,  D  =  det[P]  =  X2  -  x^. 


Consequently,  from  relations  (2.1.10a),  the  interpolation  functions  are 


<  N(x)  >  =  <  N,  (x)  N2  (x)  >  =  - - -  <  1  X  > 

(X2  -  Xi) 


X2  -Xj 
-1  1 


giving  specifically. 


X2  -  X 


Ni(x)  =  - -  N2(x)  = 


X2  -  Xi 


X  -  Xi 


X2  -  Xi 


One  verifies  explicitly  relation  (2.1.6)  in  this  case. 


The  function  f  (x)  can  now  be  approximated  with  f(x)  in  the  interval  Xj  <  x 


f(x)  =  ^  [  (X2  -  xi)fi  +  (x  -  Xi)f2 

(X2  -  x)  ^ 


To  be  specific,  choose  for  the  nodal  points: 


x,  =  0,  X,  =  _, 


(2.1.12) 

(2.1.13) 

(2.1.14a) 

(2.2.14b) 

X2: 

(2.1.15) 

(2.1.16a) 
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with  which  the  following  interpolation  functions  are  associated: 


Nj  (x)  =  -  xl  NjCx)  =  lx. 

%  2  :: 


(2.1.16b) 


Their  plots  are  shown  in  Figure  2. 
Nl,  N2 


Figure  2.  The  interpolation  fiinctions  and  Mo  associatea  with  the  two  nodal  points  Xj  =  0  and  = 
Jc/2,  respectively,  for  the  1-D  two-node  element. 


Next,  let  us  see  how  f(x)  approximates  fgvCx)  =  cos  x  with  these  two  nodes,  relation  (2.1.16b).  We  have 


fex(Xl  =  0)  =  COSXj  =  1  =  fi,  f^x  Xj  =  _  =  COSX2  =  0  =  f2. 


(2.1.16c) 


giving. 


2  (jc  1 

f(x)  = _ -  X  ,  Xj  <  X  <  X2. 

71  2 


(2.1.16d) 
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From  the  plots  (Figures  3  and  4)  we  see  that  the  agreement  is  not  perfect. 


Figure  3. 


f  [X] 


Approximation  of  f(x~)  =  cos  x  bv  two  interoolational  functions  associated  with  the  1-D  two- 
node  element.  Comparison  with  Fisure  4  shows  that  the  agreement  is  not  perfect. 


The  reason  for  this,  of  course,  is  that  the  number  of  nodal  points  is  too  small. 

To  improve  on  the  approximation,  let  us  take  three  nodes,  Xj,  X2,  and  X3.  The  polynomial  basis  is 
now  given  by  the  three-component  vector: 


<P(x)>  =  <  1,  X,  x^>. 

Following  the  previous  example  (see  also  (2.1.8a)),  one  also  writes  down  the  matrix 


[P] 


,  2 

1  Xi  Xj 

,  2 
1  X2  X2 

,  2 
1  X3  X3 


(2.1.17) 


(2.1.18) 


Unfortunately,  the  expression  for  [Q]  is  too  long  to  be  given  in  a  general  form.  However,  for  the  nodal 
points 


the  expression  is 


[Q] 


1  0  0 

__6  _8 
TC  K  K 

A.  -21  A 


(2.1.19a) 


(2.1.19b) 


From 

<N(x)  >  =  <Ni(x),  N2(x),  N3(x)>  =  <P(x)>  [Q],  (2.1.19c) 
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(2.1.19f) 


The  plot  of  f(x)  is  shown  in  Figure  6.  The  agreement  with  the  plot  for  cos  x  from  Figure  4  is  now  much 
better. 


f  [X] 


three-node  element.  The  agreement  with  Figure  4  is  now  excellent 


Now  we  take  an  example  in  the  two-dimensional  (2-D)  space.  In  order  to  simplify  the  discussion,  we 
specify  the  nodal  points  from  the  beginning: 

Tl  =  (1,1),  Tj  =  (-1,1),  T3  =  (-1,-1),  ?4  =  (1,-1),  (2.1.20) 

where  in  general  f  =  (x,y).  The  element  associated  with  these  four  nodal  points  is  referred  to  as  a  2-D 
four-nodal  quadrilateral  element,  m  =  4.  Its  plot  is  exhibited  in  Figure  7. 
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(Because  of  the  way  the  nodal  points  are  chosen,  this  and  other  elements  like  this  one  are  also  referred 
to  as  "reference  elements").  As  far  as  the  polynomial  basis  functions  are  concerned,  they  are  defined 
through  the  four-component  vector  (n  =  4): 

<P(r)>  =  <  1,  X,  y,  xy>.  (2.1.21) 

Following  the  general  prescription  for  constructing  the  interpolation  functions,  with  the  help  of 
<P(r)>  =  <1,  X,  y,  xy>  evaluated  at  positions  of  nodal  points,  we  first  write  down  the  four-by-four 
matrix 


[P] 


1111 
1-1  1-1 
1-1-1  1 
1-1  1-1 


(2.1.22a) 
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which  when  inverted  yields 


[Q] 


1  1  1 

-1  -1  1 

1  -1  -1 

-1  1  -1 


(2.1.22b) 


Now  applying  now  relations  (2.1.22b)  and  (2.1.21)  to  relation  (2.1.10a)  one  obtains  the  following 
expression  for  the  interpolation  functions: 


N,  (x,y)  =  i.  (1  +  X  +  y  +xy),  N2(x,y)  =  .i(l-x+y  -  xy), 

'4  4 

N3(x,y)  =  .i(l-x-y  +  xy),  N4(x,y)  =  +  x-  y-  xy).  (2.1.22c) 

4  4 


In  Figure  8,  the  first  of  these  four  functions  is  displayed.  The  other  three  are  obtained  by  rotating  the  x-y 
plane  about  z-axis  clockwise  by  2^,  and  3^  angles,  respectively. 


To  demonstrate  how  things  woik  in  the  2-D  space,  we  choose  for  the  exact  function  the  expression; 

fex(^>y)  =  exp[-(x  +  y)].  (2.1.23a) 

From  relations  (2.1.20)  the  nodal  parameters  are: 

=  e~2  =  0.13534,  f2  =  1,  fj  =  e^  =  7.38906,  f4  =  1;  (2.1.23b) 

yielding  on  this  element  for  exp  [-(x  +  y)]  the  appropriate  expression: 

f(x,y)  =  0.135  Ni(x,y)  +N2(x,y)  +  7.389  N3(x,y)  +  N4(x,y).  (2.1.23c) 

The  plots  of  and  f  are  exhibited  on  Figures  9  and  10,  respectively. 
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Figure  10.  Approximation  of  exp  r-(x  +  v)1  by  the  four  interpolation  fiinctions  associated  with  the 
quadrilateral  element  from  Figure  8. 

The  agreement  is  quite  good,  taking  into  account  that  we  used  only  four  nodal  points,  which  is  not  a  large 
number  for  the  2-D  space.  If  the  number  of  nodal  points  is  increased  at  least  by  one,  the  agreement  would 
be  perfect. 

Of  course,  there  are  many  more  elements  that  one  can  consider,  not  just  in  1-  and  2-D  spaces,  but, 
in  fact,  in  an  arbitrary  dimensional  space.  However,  such  a  study  is  beyond  the  scope  of  this  repoit 

Finally,  we  illustrate  briefly  how  to  solve  a  partial  differential  equation  by  finite  element  method.  To 
be  specific,  we  treat  Poisson’s  equation  in  the  2-D  space  defined  over  a  square  region  (which  is  actually 
the  quadrilateral  element  from  the  previous  example)  with  a  constant  surface  "charge  density"  a  (consult 
Figure  11  for  details). 
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y 


Figure  11.  Square  region  with  a  constant  surface  charge  a  generating  the  solution  f  from 
equation  in  2-D  space.  At  the  boundary,  the  solution  f  is  required  to  vanish. 

We  write  the  differential  equation  as  follows: 

2(0  *  o  =  .0-0, 

dx^  dy^ 


where 


a  0:  -1  <  X  <  1,  -1  <  y  <  1;  a  =  0  elsewhere. 

One  easily  sees  fiom  equation  (2.1.24a)  that  the  solution  possesses  the  following  symmetries: 

f(x,y)  =  f(y,x);  f(x,y)  =  f(-x,y)  =  f(-x,-y)  =  f(x,-y). 
Furthermore,  we  also  require  that  f  satisfy  these  boundary  value  conditions: 


Poisson 


(2.1.24a) 


(2.1.24b) 


(2.1.24c) 
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f(±l,y)  =  f(x,±l)  =  0. 


(2.1.24d) 


Because  of  the  symmetry  and  the  boundary  value  conditions,  representation  of  f  in  terms  of  a  just-derived 
four-node  quadrilateral  element  would  not  be  enough;  one  would  have  to  constract  at  least  a  five-node 
element  However,  rather  than  do  that,  let  us  actually  exploit  the  symmetry  of  the  problem,  and  instead 
of  interpolation  functions,  use  the  polynomial  functions  as  discussed  at  the  beginning  of  this  section.  To 
simplify  the  discussion,  we  use  just  two  of  them;  the  simplest  ones  that  satisfy  conditions  of  symmetry 
and  the  boundary  values  are: 

Pi(x,y)  =  (x^  -  1)  (y^  -  1);  P2(x,y)  =  (x^  +  y^)  Pi(x,y).  (2.1.25a) 


Their  plots  are  shown  on  Figures  12  and  13. 


Figure  12.  Plot  of  the  lowest  order  polynomial  function  with  the  same  symmetries  as  the  function  f  from 
Figure  11. 
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Figure  13.  Plot  of  the  second  lowest  order  polynomial  fiinction  with  the  same  symmetries  as  the 
function  f  from  Figure  11. 

These  functions  allow  us  to  express  the  f  in  terms  of  two  generalized  parameters  of  approximation  a^  and 

ai: 


f(x,y)  =  <  Pi(x,y)  P2(x,y)  > 


Pi(x,y)  aj  +  P2(x,y)  a2. 


(2.1.25b) 


One  notices  that  now  the  differential  equation  (2.1.24a)  decomposes  as 

££(0  =  ai££(Pi)  +  a2££(P2),  9e(Pi)  =  2(x2  +  y^  _  2), 
c£(P2)  =  2(6x2  _  1)  (y2  _  1)  ^  2(6y2  -  1)  (x^  -  1) 

+  2(x'^  -  x2)  +  2(y^  -  x2).  (2.1.25c) 
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The  problem  now,  of  course,  is  to  find  and  32- 

In  order  to  determine  and  82,  we  choose  two  collocation  points  (points  inside  the  square  where 
o;6  0): 


Ti  =  (0,0),  T 


(2.1.26a) 


and  define,  in  general,  the  weight  functions 

Wi  =  (g?(f)  +  0)y,  =  (aj^OPi)  +  a2  9e(P2)  +  o)y.;  Wj  =  0.  (2.1.26b,c) 


The  statement  (2.1.26c)  means  that  the  differential  equation  is  satisfied  exactly  at  the  collocation  points. 
One  obtains  explicitly 


Wi  =  -4  at  +  4  82  +  o  =  0,  W2 


+  a  =  0, 


(2.1.27a) 


with  solutions 


aj  =  0.2976a,  82  =  0.0476a.  (2.1.27b) 

This  gives 

f(x,y)  =  0.2976a  (x^  -  1)  (y^  -  1)  +  0.047a  (x^  +  y^)  (x^  -  1)  (y^  -  1).  (2.1.27c) 
The  function  f  is  shown  in  Figure  14  for  a  =  1. 
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Our  result  compares  very  well  with  the  one  obtained  from  the  Fourier  series  expansion.  Specifically, 
from  there  one  obtains  for  f  at  T  =  0  the  value: 

Fourier  Series:  f(f  =  0)  =  0.2976o.  (2.1.27d) 

On  the  other  hand,  directly  from  solution  (2.1.27c),  we  have 

f(T  =  0)  =  ai  =  0.2976a,  (2.1.27e) 

which  is  practically  the  same  as  the  Fourier  series  result. 

Although  rather  simple,  the  few  examples  chosen  here  illustrate  clearly  the  power  of  the  finite  element 
method.  The  interested  reader  can  find  more  involved  examples  in  a  rather  extensive  survey  of  the 
literature,  some  of  which  is  listed  at  the  end. 
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2.2  Finite  Difference  Method.  In  the  finite  difference  method,  as  the  name  would  indicate,  one  uses 


differences  of  functions,  variables,  and  the  like,  in  order  to  obtain  approximate  solutions  to  mostly 
differential  equations.  In  finding  munerical  solutions,  the  difference  operators  are  used  in  a  manner  similar 
to  that  of  differential  operators  in  the  differential  equations  (Hovanessian  1976). 

In  order  to  define  various  difference  operators,  we  consider  a  function  y  =  f(x)  as  shown  in  Figure  15 
with  equal  intervals  Ax  for  an  independent  variable  x. 


Figure  15.  Definition  of  the  equal  interval  Ax  from  which  the  forward,  the  backward,  and  the  central 
differences  of  v(x)  are  defined. 

There  are  three  kind  of  differences:  the  forward,  the  backward,  and  the  central  difference,  respectively. 
The  first  forward  difference  of  yj  is  defined  as 

^yi  =  yi  + 1  -  yi>  (2.2.1) 
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where  fomally  A  is  the  forward  difference  operator.  Qeaily,  for  the  second  forward  difference  operator 
we  have 


A^yj  =  ACAyj)  =  Ayj  ^  1  -  Ayj  =  yj  ^  2  "  2yi  +  i  +  yi-  (2.2.2) 

By  similar  reasoning  one  obtains  for  the  third  forward  difference  operator  expression 

A^y;  =  ACA^yj)  =  yi  +  3  "  3yi  +  2  +  ^yj  ^  1  -  yj.  (2.2.3) 

These  two  equations  clearly  exhibit  the  operational  nature  of  A. 

The  first  backward  difference  of  yj  is  defined  as 

Vyj  =  yj  -  yi_i,  (2.2.4) 

where  formally  V  is  the  backward  difference  operator.  The  second  backward  difference  operator  can  be 
obtained  from  the  first: 

V^yj  =  V(Vyi)  =  Vyj  -  Vyj  _  ^  =  yi-2yi  _  1  +  yi  _  2-  (2-2.5) 

Similarly,  the  third  backward  difference  becomes 

V^yj  =  V(V2yj)  =  yi  -  3yi  _  1  +  3yj  _  2  -  yi  _  3.  (2.2.6) 

The  first  central  difference  of  yj  is  defined  as 

5yi  =  y.  1  -  y.  1 ,  (2.2.7) 

‘"t  t 
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where  5  formally  denotes  the  central  difference  operator.  The  second  central  difference  operator  is 
obtained  from  the  first: 


5^yi  =  6(6yi)  =  5^’ 


V  2  2^ 


=  yi  + 1  -  2yi  +  yi  _  1- 


y.  1  1  -  y.  1  1 

^  *  -K-  *  ^  *  -ir  ~  -TT 

T  T  T  T 


y.  1  1  -  y.  1  1 

^  "  T  ^  T  *  "  T 


Similarly,  third  and  fourth  central  difference  operators  are  obtained: 


5^yi  =  6(6^yi)  =  yj  ^  2  "  + 1  -  I  ~  -  1’ 

2  2  2  2 

6^yi  =  6(6^yi)  =  yi  +  2  “  ^yi  +  2  +  ^yj  -  4yi  _  ^  +  yi  _  2- 


Denoting  generically  with  the  three  difference  operators,  one  verifies  directly  that 


f22(Qy.)  =  n(n2y.). 


By  induction  one  further  obtains 


=  Q"  ^  =  1, 


and  very  important  distribution  laws 


Q(yj  +  Zj)  =  f2yi  +  Qzj, 


Qconstyj  =  constfJyj. 


(2.2.8) 

(2.2.9) 

(2.2.10) 

(2.2.11) 

(2.2.12) 

(2.2.13) 

(2.2.14) 
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Next,  we  wish  to  establish  the  relationship  between  difference  and  differential  operators.  To  this  end, 
define  first  the  symbolic  difference  operator  E  as  follows: 


Eyi  =  yi*i. 


E"^ 


yi 


(2.2.15) 


One  should  notice  that  these  relations  are  consistent  with  each  other.  Some  important  relations  involving 
E  are: 


E  const  =  const  E, 


EO  =  O  E,  E  "O  =  QE  ",  (2.2.16) 

where  const  and  n  are  independent  of  x.  Using  these  relations,  one  obtains  difference  operators  A,  V,  and 
5  in  terms  of  E 


Ayj  =  (E  -  l)yj  >-*•  A  =  E  -  1, 


Vyj  =  (1  -  E  ^)yj  V  =  1  -  E 


1  _  1 
5  =  e"^  -  E 


(2.2.17) 


The  relation  between  the  forward  difference  operators  and  the  differential  operator  is  obtained  by 
writing  Taylor’s  series  expansion  of  the  function  f(x  +  h)  about  x. 


f(x  +  h)  = 


E 


n  =  0 


ii^D  "f(x)  =  e*‘°f(x),  D  =  -1. 

n!  dx 


(2.2.18) 
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But 


Ef(x)  =  f(x  +  h)  »-»  E 


(2.2.19) 


Now  we  are  ready  to  find  a  relationship  between  difference  operators  and  D.  Starting  with  A,  from 
relation  (2.2.17)  we  have 


A  =  E-  I=  e*’®-1,  D  =  .iln(l  +  A). 


(2.2.20) 


The  useful  relations  are  obtained  when  one  expands  (2.2.20),  yielding: 


A  =  hD  +  —D^  +  +  ... 


2! 


3! 


(2.2.21) 


and 


D  =  _L 
h 


^  _  A^  ^  A^  _  A^  ^ 

T  ^  ~  ~T  ^ 


(2.2.22) 


Higher-order  derivative  and  forward  difference  operators  are  obtained  by  taking  various  powers  of  (2.2.20), 
resulting  in  the  following  second-  and  higher-order  operator  equations: 


A^  =  h^D^  +  h^D^  +  -Lh'^D^  +  ..., 

12 


A^  =  h^D^  +  Ih^D^  +  £h^D^  + 

2  4 


(2.2.23a) 

(2.2.23b) 


=  J-fA^  -  A^  +  Ha^  -  1a^  +  ... 

12  6 


h\ 


(2.2.24a) 
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(2.2.24b) 


=  J_ 


A’  -  1a*  -  Ia’  * 
2  6 


h" 


V  -  2^^  +  JJ-A”  -  ... 


17  a6 


(2.2.24c) 


Equations  (2.2.21)  to  (2.2.24)  can  be  used  to  find  the  differentials  of  functions  in  tenns  of  available 
functional  values.  For  example,  neglecting  errors  of  order  h  and  higher,  D^yj  can  be  obtained  firom 
(2.2.24a).  Dropping  A^  and  higher-order  terms  from  this  equation  and  using  (2.2.2),  one  obtains 


D^yi  =  .lA^y,  -  ^(yi  .  2  -  2yi  .  j  *  yj. 


(2.2.25a) 


A  faster  way  to  obtain  this  result  is  to  use  (2.2.20): 

D^y,  =  ’  A^y,  =  ’  (1  -  E)2y,  -  4  (^i  ’  *  1  Xi  •  2)- 

\r  h^  h^ 

The  backward  difference  operator  can  be  obtained  in  terms  of  the  differential  operator  D  with  the  aid 
of  equation  (2.2.17).  In  fact,  from  the  relations  (2.2.17)  and  (2.2.19)  one  has 

V  =  1  -  E-l  =  1  -  e-‘'°,  D  =  -lln(l  -  V).  (2.2.26,27) 

h 


Taylor’s  expansion  yields  for  V  the  result 


V 


(2.2.28a) 


Taking  higher  powers  of  relations  (2.2.21a,b)  we  obtained  higher-order  difference  and  derivative  operators. 
These  are  summarized  as  follows: 
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(2.2.28b) 


-..., 

12 

-  Ih^D^  +  Ih^D^  -  ...;  (2.2.28c) 

2  4 


f 


(2.2.29a) 


+ 


+  iiv'*  +  £v^  + 
12  6 


{2229b) 


= 


+  +  Zv^  + 

2  4 


{2229c) 


+  2V^ 


(2.2.29d) 


The  equations  relating  the  central  difference  operator  5  and  differential  operator  D  are  obtained  using 
equations  (2.2.15)  and  (2.2.19). 


1  _i 

6  =  -  E  Z  E  =  e  d  =  llnE. 

h 


(2.2.30) 


In  addition,  it  is  customary  to  define  the  mean  operator  p  as  follows: 


P 


1 

•  T 


+  E 


0 

2 


(2.2.31) 
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This  implies 


p5  =  1(e  -  E-i)  =  e-^^)  =  sinh(hD). 

2  2 


Applying  the  difference  operator  p5  to  y^,  we  obtain 

0i5)yi  =  l(yi  ^  1  -  yi  _  i). 


From  (2.2.31),  we  have 


The  Taylor  series  expansion  of  the  right-hand  side  of  (2.2.32)  yields 


p6 


hD  + 


h^D^ 


h^D^ 

“W 


The  general  expression  for  5",  where  n  =  1,  2,  3, ...,  can  be  obtained  from  relation  (2.2.30). 


hD  _hD 

6  =  e~^  -  e  =  2  sinh  ilB., 

2 

6"  =  2"  sinh" 

2 


By  similar  reasoning,  we  obtain  from  (2.2.31)  for  n  =  1, 2,  3, ..., 


(2.2.32) 


(2.2.33) 


(2.2.34) 


(2.2.35) 


(2.2.36) 
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(2.2.37) 


P 


=  cosh 


hD 


«  n  hD 
=  cosh - 


Using  the  combinations  of  the  aforementioned  relations,  we  can  obtain  expressions  for  higher  powers 
of  central  difference  operators,  as  for  example 


p5^  =  (ji5)  5^  =  sinh  hD  ir  sinh 


2  hD 


h^D^ 


h^D^ 

40 


Similarly,  from  (2.2.36)  one  obtains 


5-4  .4rx4  h^D^  h^D* 

5  =  h  D^  + _  + _  + 

6  80 


(2.2.38) 


(2.2.39) 


The  differential  operators  D  are  obtainable  in  terms  of  central  difference  operators  6.  For  example, 
from  (2.2.32)  and  (2.2.35)  we  have 


r 


D  =  JL  sinh"^u5  =  -L 
h  h 


^  6  30 


(2.2.40) 


where  the  right  side  is  obtained  from  the  Taylor  series  expansion  of  sinh  Similarly, 


D^  =  JL 


12  90 


(2.2.41) 


d3  =  JL 


V  4  120  , 


(2.2.42) 
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(2.2.43) 


6 


) 


Although  not  exclusively  true  for  all  cases,  when  numerically  solving  differential  equations,  one  tends  to 
express  differential  operators  D  in  terms  of  difference  operators  5. 

As  an  example  of  application  of  the  aforementioned  equations  with  difference  operators  6,  we  seek 
the  numerical  solution  of  the  differential  equation 

2 

£2  +  +  y  +  0;  y(t  =  0)  =  1,  (dy/dt)  (t  =  0)  =  0.  (2.2.44a,b) 

dt2  dt 

Using  differential  operator  D  =  d/dt,  this  differential  equation  translates  into 

D^yj  +  Dyj  +  y^  =  0;  yo  =  1,  Dyg  =  0.  (2.2.45a,b,c) 

Denoting  At  =  h  and  containing  relations  (2.2.30)  and  (2.2.31)  with  relations  (2.2.40)  to  (2.2.43),  we  obtain 

Dy  =  ^  ~  D^y.  =  ^  ~  ^  e  =  0(h^),  (2.2.46a,b) 

2h 


where  e  stands  for  "error"  and  denotes  the  order  in  h  of  terms  that  are  neglected.  Substituting  relations 
(2.2.46)  into  (2.2.45a),  we  obtain  a  recursive  relation  for  yj  ^  j,  as  follows: 


yi  +  1  = 


1 


2  +  hL 


(4  -  2h2)  yj  +  (h  -  2)  yj.  1 


(2.2.47) 


This  equation  gives  yj  + 1,  providing  the  values  of  yj  and  yj.j  are  available:  for  these  the  initial  conditions 
will  be  helpful.  The  condition  (2.2.45c)  gives 
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yi  “  y-i 

Dy«  - -4ipi  '  o;  y_.  =  y,  . 


(2.2.48) 


We  inserted  y_i  in  order  to  make  the  derivative  at  yo  symmetrical.  Formally,  y_i  can  be  related  to  the 
second  constant  of  integration.  In  any  case,  setting  i  =  0  in  (2.2.47)  and  using  yg  =  1  and  y_i  =  yj,  we 
obtain 


yj  =  .^[(4  -  2h^)  (h  -  2)yi 
2  +  h*- 


(2.2.49a,b) 


Selecting  the  value  for  the  increment  h,  we  obtain  a  numerical  value  for  yj.  Thus  we  will  have  yg  and 
yj,  and  we  can  proceed  with  the  calculation  of  y2  by  setting  i  =  2  in  (2.2.47).  Similarly,  we  can  calculate 
y3,  y^,  and  so  on.  The  values  yg,  y^  y>i,  y-^,  ...  will  represent  the  solution  of  the  differential  equation 
(2.2.44)  at  discrete  time  intervals  of  i  At  =  ih  where  i  =  0,  1,  2,  3 . 

Already  we  see  the  possible  problems  with  the  finite  difference  method  solutions  of  differential 
equations.  If  At  is  too  small,  a  great  number  of  iterations  will  be  required  to  obtain  the  value  of  y  at  a 
given  time.  On  the  other  hand,  if  At  is  too  large,  since  the  error  is  of  the  order  of  (At)^  =  h^,  the  error 
in  each  application  of  (2.2.47)  may  accumulate,  resulting  in  large  errors  in  the  solution. 

Next  we  discuss  a  simple  case  of  numerical  solution  of  a  differential  equation  starting  from  the  Taylor 
series  expansion.  Consider  the  first-order  differential  equation 

y'  =  -^  =  f(x,y),  (2.2.50) 


with  the  initial  condition  y  =  yg  at  x  =  Xg.  The  Taylor  series  expansion  of  y  about  Xg,  using  an  increment 
Ax,  yields 


y(xo  +  Ax)  =  yg  +  Axyg  +  ^(Ax)2yo  + 


2.2.51 
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Denoting  Xj  =  Xq  +  Ax  and  Xj  =  Xj  +  Ax  with  corresponding  yj  =  y(xo  +  Ax)  and  y2  =  yCx^  +  Ax),  we 
may  write 


Yi  =  Yo  +  Axf(xo,yo). 
y2  =  Yi  +  Axf(xi,yi), 


yn  +  1  =  yn  +  Axf(xn«yn)- 


(2.2.53) 
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The  error  is  simply  the  next  term  in  the  expansion  (2.2.51): 


.  _  (Ax)\,/  _  (Axf 
2  ^  2 

The  error  in  this  case  is  directly  proportional  to  the  square  of  the  increment  Ax  and  to  the  combination 
of  function  f(x,y)  and  its  partial  derivatives,  with  respect  to  x  and  y  at  the  position  (x,y)  =  (Xj.yj),  i  =  1, 
2,  3, ...,  a 

Example:  Consider  the  numerical  solution  of  the  first-order  differential  equation 

=  y^  =  -y  with  y  =  1  at  x  =  0.  (2.2.55a) 

dx 


df 

dx 


3y 


(2.2.54) 


The  solution  is  given  exactly  as 


y  =  e"\  (2.2.55b) 

Denoting  with  Ax  =  h,  from  relations  (2.2.53),  one  immediately  sees  that 


y„  =  1  -  nh. 


(2.2.56) 


which  is  a  good  approximation  of  exact  solution  (2.2.55b)  if  x„  s  nh  <  x 


2 

n  * 


Another  way  to  solve  this  differential  equation  numerically  is  to  use  the  finite  difference 
representation.  Namely,  the  derivative  of  y  at  Xj,  is  equal  to  the  slope  of  the  curve  of  Figure  16  at  this 
point.  This  slope  can  be  represented  by  either 

yn  - 
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or 


(2.2.58) 


yn  =  (yn.l-yn-l)^^. 


where,  of  course,  h  =  Ax.  One  recognizes  equations  (2.2.57)  and  (2.2.58)  as  the  forward  and  central 
difference  equations  of  the  first  derivative.  Therefore,  consistent  with  relation  (2.2.55a)  one  may  write 


(yn  *  1  -  yn)/h  =  -yn^  yo  =  1-  h  =  ^, 


(2.2.59) 


(yn  >  1  -  yn-i)/2h  =  -y„;  yo  =  1.  yi  =  e  h  =  1. 


(2.2.60) 


Formulation  (2.2.59)  is  the  forward  difference  formulation  with  increment  h  =  0.25.  Formulation 
(2.2.60),  however,  is  the  central  difference  formulation,  also  with  h  =  0.25.  Now,  in  formulation  (2.2.59), 
we  can  set  n  =  0  and  solve  for  y^.  Formulation  (2.2.60)  will  require  both  yg  and  y^  to  start  the  problem 
and  solve  for  ^2-  ™ost  cases,  y^  can  be  evaluated  by  using  the  Taylor  series  expansion. 

Of  course,  there  are  other  finite  difference  methods  to  numerically  solve  differential  equations;  one 
such  method  is  the  famous  Runge-Kutta  method.  Discussing  all  these  methods  is  beyond  the  scope  of  this 
report. 

2.3  Finite  Volume  Method.  Finally  we  address  briefly  the  finite  volume  method.  In  the  finite 
volume  method  (Hyman,  Kn^p,  and  Scovel  1992),  the  average  values  of  a  function  over  local  me^  cells 
are  taken  as  unknowns;  discrete  approximations  of  the  divergence,  gradient,  and  rotor  (curl)  operators  are 
defined  using  general  forms  of  Stoke ’s  theorem. 

To  get  a  feel  for  the  finite  volume  method,  consider  a  physically  motivated  system  of  partial 
differential  equations  derived  from  a  limiting  process  applied  to  integral  equations.  For  example,  a 
quantity  p  (charge  density  or  fluid  density)  is  conserved  under  the  flow  of  a  conservation  law  if  the 

amount  of  p  contained  in  any  fixed  volume  is  due  entirely  to  the  flux  of  j  (current  density)  across  the 
boundary  o(f2)  of  the  volume  £2;  one  should  realize  that  here  the  terms  "volume"  and  "surface"  are  not 
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necessarily  restricted  to  just  the  three-dimensional  (3-D)  space.  For  the  sake  of  simplicity,  however,  here 
we  do  deal  in  the  3-D  space,  and  the  conservation  law  can  be  expressed  in  integral  form  as 


_d 


o(£2) 


dS  , 


(2.3.1) 


where  dSis  the  surface  element  vector  whose  direction  is  outward  normal  to  the  boundary.  Moving  the 
time  derivative  under  the  integral  sign  and  applying  the  divergence  (Gauss)  theorem,  equation  (2.3.1)  can 
be  rewritten  as 


p  +  V  -jjdV  =0. 


(2.3.2) 


If  we  let  the  volume  shrink  to  a  point,  we  obtain  the  partial  differential  equation 

atP  +  V-j=0,  (2.3.3) 


providing  that  at  every  point  p  and  j  are  differentiable. 

When  numerically  solving  equation  (2.3.1)  it  is  natural  to  stop  the  limiting  process  at  the  local  mesh 
spacing  and  solve  (2.3.2)  where  the  control  volumes  are  the  local  mesh  cells.  Specifically,  in  a  small 
time-invariant  3-D  control  volume  Q,  we  rewrite  the  integral  equation  (2.3.2)  as 


(P).v3  +  (V  •  j)..3  =  0, 

where  (’’OavS  ’  )av3  represent  a  3-D  cell  average;  the  divergence  theorem  is  applied  in 

computing  (V  •  j  by  acting  on  face-average  normal  components  of  fluxes: 


35 


’DavS  •=  YqJ 


dVol 


I _  f  j  •  d5 

(0\  j  •' 


Vol  (Q) 


o(£2) 


Vol  (Q)  ? 


av2. 


(2.3.4) 


Here  the  time-independent  control  volume  Q  is  bordered  by  mesh  points.  Its  boundary  is  the  union  of 
I  (i  =  1, 2, ....  I)  distinct  pieces,  a(f2)  =  UaCfij).  Furthermore,  Aj  is  the  area  vector  of  the  i  th  piece,  and 

(ji)av2  niay  be  interpreted  as  the  normal  component  of  j  over  the  i  th  piece  given  by 


ai)av2  =4-  /  j  (2-3.5) 


The  central  idea  behind  the  finite  volume  method  is  to  accurately  approximate  (ji)av2  use  relation 
(2.3.4)  to  define  a  discrete  approximation  of  (V  •  j}av3,  which,  in  turn,  defines  8(p)ay3/dt. 

In  fact,  the  finite  volume  method  can  be  directly  applied  to  time;  integrating  (2.3.2)  from  to  j 
we  obtain 

•n  ♦  1 

J  dt  J  (9tp  +  V  •  j)  dVol  =  0,  (2.3.6) 


or,  using  (2.3.4),  one  further  obtains. 
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(pt;‘  =  (pt3-  j  (v-j)„3di 

tn 

»  1 

■  <*’  "  Wffi)  ?  /  <2.3.7) 

The  fluxes  associated  with  (ji)av2  assumed  to  be  known  as  a  function  of  time;  they  might  be 
approximated  directly  by  incorporating  past  time  levels  in  standard  linear  multistep  methods  (Hyman, 
Knapp,  and  Scovel  1992);  of  course,  there  is  nothing  wrong  in  evaluating  these  fluxes  by  the  finite 
difference  of  the  finite  element  method  if,  for  example,  j  is  given  by  some  other  equation  (the 
Schroedinger  equation,  for  instance). 

This  was  a  simple  "introductory"  example  intended  to  give  a  "flavor"  and  the  basic  idea  of  the  finite 
volume  method.  It  is  beyond  the  scope  of  this  report  to  go  into  more  details  of  this  rather  new  and 
complex  method,  which,  at  least  presently,  can  be  found  only  in  research  literature. 

3.  COMPARISONS  OF  METHODS 

Today,  both  finite  element  and  finite  difference  methods  are  used  rather  extensively  in  numerically 
solving  the  most  complex  problems,  as  in  calculating  the  electromagnetic  fields  for  complex  structures  or 
in  calculating  the  electromagnetic  field  effects  on  structures  with  very  complex  geometries.  However,  it 
is  the  finite  element  method  that  is  today  considered  to  be  a  superior  method,  particularly  when  the 
question  of  efficiency  is  posed  for  solving  a  variety  of  practical  problems.  One  of  the  main  attractions 
of  this  method  is  that  once  the  method  is  set  up,  that  is,  the  computer  program  is  written  which  generates 
the  interpolation  functions,  it  can  be  used  for  other  problems  by  changing  the  input  data.  On  a  historical 
note,  the  method  was  first  developed  in  1956  for  the  analysis  of  aircraft  structural  problems.  Today,  so 
much  work  is  done  with  the  finite  element  method  that,  in  fact,  this  method  has  become  one  of  the 
research  areas  of  applied  mathematics  itself 

The  finite  difference  method,  of  course,  is  much  older  than  the  finite  volume  method.  In  fact,  it  is 
as  old  as  differential  calculus  itself  This  means  that  it  was  in  use  well  before  computers  became 
available.  However,  each  advance  in  computer  technology  makes  the  finite  difference  method  more 


37 


interesting  in  that  more  and  more  complex  problems  can  be  treated.  In  fact,  being  "simpler"  than  the  finite 
element  method,  the  finite  difference  method  is  quite  often  the  preferred  method  when  the  accuracy  of 
results  is  not  essential. 

In  fact,  rather  recently,  at  the  Stanford  Linear  Accelerated  Center,  where  the  solutions  of  Maxwell’s 
equations  for  electromagnetic  fields  propagating  in  cavities  with  symmetric  periodic  structures  and 
structures  with  axial  symmetry  were  sought,  it  was  foimd  that  the  finite  element  method  is  much  more 
accurate  than  the  finite  difference  method  (Nelson  1992, 1993).  Denoting  genericaUy  with  F  the  solution 
for  an  electromagnetic  field,  one  finds  that  the  smallest  values  for  |5F/F|  are 

|5F/F|  =  10“^  10"^ 

for  the  finite  difference  and  finite  element  methods,  respectively.  Here  6F  represents  the  deviation  from 
the  experimentally  measured  field  F.  Furthermore,  taking  into  account  that  the  finite  element  method  is 
capable  of  finding  rather  accurately  the  solutions  for  electromagnetic  fields,  even  when  the  structures  have 
sharp  edges,  one  sees  that,  overall,  the  finite  element  method  is  superior  to  the  finite  difference  method. 
However,  the  finite  element  method,  because  it  follows  rigorous  procedures,  can  be  rather  complex  when 
deriving  numerical  algorithms  for  solutions,  particularly  when  compared  to  the  finite  difference  method. 
Hence,  when  analyzing  technical  and  scientific  problems,  these  two  methods  should  not  be  considered  as 
mutually  exclusive,  but  rather  both  of  them  should  be  used,  depending  on  needed  accuracy  and  urgency 
of  solutions. 

As  far  as  the  finite  volume  method  is  concerned,  its  accuracy  depends  on  specifics  of  other 
computational  methods  that  are  used  after  the  generalized  Stoke’s  theorem  has  been  applied.  The  thing 
to  emphasize  here  is  that  the  generalized  Stoke’s  theorem  in  the  finite  volume  method  effectively  reduces 
the  number  of  computational  dimensions;  a  solution  that  is  being  sought  at  a  small  volume  (which  mimics 
a  point)  is  reduced  with  the  help  of  the  generalized  Stoke’s  theorem  to  computations  over  sides  bordering 
a  small  volume. 
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4.  POTENTIAL  APPLICATIONS,  DISCUSSIONS,  AND  CONCLUSIONS 


Finite  volume  and  finite  difference  methods,  no  doubt,  can  be  used  to  obtain  adequate  approximate 
solutions  of  electromagnetic  fields.  In  particular,  one  should  be  able  to  predict  the  effects  of  nuclear-burst¬ 
generated  electromagnetic  effects  on  defense  systems,  such  as  the  scattering  of  the  electromagnetic  pulse 
off  combat  vehicles  as  a  whole  and  electronic  components,  sensors,  etc.,  in  particular. 

One  of  the  newer  subjects  where,  in  particular,  the  finite  element  method  could  be  used  is  the 
statistical  modeling  of  the  Gulf  War  Syndrome.  The  aim  here  is  to  obtain  the  probabalisitic  formulation 
of  the  syndrome  itself.  Specifically,  for  each  unit  associated  with  a  particular  geographic  location  based 
on  data  from  the  U.S.  Army  &  Joint  Services  Support  Group  and  the  Environmental  Epidemiological 
Service  from  the  Veterans  Administration,  the  soldier  probability  distribution  function  of  symptoms  (here 
symptom  means  symptom  and  illness)  would  be  constructed  with  the  help  of  six  finite  element 
interpolation  functions  (corresponding  to  six  discrete  symptoms).  Because  of  these  interpolation  functions, 
the  soldier  probability  distribution  function  becomes  an  analytical  function  in  symptom-variable.  This, 
in  turn,  allows  one  to  calculate  the  average  soldier  S)miptom  as  a  function  of  the  continuous  symptom- 
variable,  allowing  us  to  find  the  range  of  symptoms  that  very  ill  soldiers  have.  Turning  now  to  the 
database  we  can  rank  specific  symptoms  according  to  the  frequency  of  their  occurrence.  The  six  most 
frequent  ones  define  the  syndrome  which  now  the  medical  doctors  can  study  to  determine  the  common 
causes  of  the  Gulf  War  Syndrome  for  each  particular  geographic  region. 

Another  new  project  where  actually  the  finite  volume  method  could  be  employed  is  the  "Rotor  Wash 
Model"  project.  Here  the  idea  is  to  develop  analytical  techniques  which  would  replicate  rotorcraft 
aerodynamics  at  specified  positions  around  the  helicopter  fuselage  in  order  to  predict  the  infiltration  of 
diverse  kinds  of  agents  on  and  into  helicopters.  It  appears  natural  to  model  the  flow  of  agents  with  the 
finite  volume  method  by  employing  the  Stoke ’s  theorem  when  studying  the  flow  of  agents. 

Finally,  we  notice  that  "Aerosol  and  Vapor  Infiltration  Modeling"  is  another  important  area  where, 
in  addition  to  the  finite  element  and  finite  difference  methods,  the  finite  volume  method  could  play  a  very 
important  role.  Here,  one  is  striving  to  develop  a  predictive  methodology  that  could  determine  levels  of 
vapor  concentration  and  liquid  deposition  in  and  on  a  variety  of  stmctures.  The  fluid  dynamics,  with 
which  one  is  dealing  here,  are  natural  for  utilization  of  the  finite  volume  method.  This  again  comes  from 
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the  fact  that  Stoke’s  theorem,  on  which  the  finite  volume  method  is  based,  enters  very  naturally  into  the 
importance  of  the  description  of  fluid  motion. 

Hopefully,  the  few  examples  mentioned  here  show  the  importance  of  the  methods  briefly  exposed  in 
the  main  body  of  the  text.  It  appears  that  the  finite  volume  method  is  not  yet  widely  accepted  by  the 
computing  community.  The  reason  for  this  is  that  the  finite  volume  method,  when  compared  to  the  other 
two  methods,  is  not  a  "pure"  method  but  rather  a  "composite"  method  consisting  of  Stoke’s  theorem  and 
the  finite  element  and/or  the  finite  difference  methods.  Judging  by  the  amount  of  literature  available  on 
various  methods  (see  the  Appendix  for  details)  it  is  evident  that  the  finite  element  method  is  far  the  most 
popular  and  perhaps,  the  most  important  approximation  method  for  evaluations  in  vulnerability/  lethality 
analysis. 
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DR  WADE 

WSMR  NM  88002-5513 

1  ARMY  RESEARCH  LABORATORY 

ATTN  AMSRL  SL  E 
MR  MARES 
WSMR  NM  88002-5513 

1  ARMY  TRADOC  ANL  CTR 

ATTN  ATRC  W 
MRKEINTZ 
WSMR  NM  88002-5502 

1  ARMY  TRNG  &  DOCTRINE  CMND 

ATTN  ATCD  B 
FT  MONROE  VA  23651 

ABERDEEN  PROVING  GROUND 

1  CDR  USATECOM 
ATTN;  AMSTE-TA 

2  CDR  USAMSAA 
ATTN:  AMXSY-ST 

AMXSY-D 

4  DIR  USARL 

ATTN:  AMSRL-SL,  J  WADE  (433) 

AMSRL-SL-I  M  STARKS  (433) 
AMSRL-SL-C,  W  HUGHES  (E3331) 
AMSRL-SL-B,  P  DEHZ  (328) 

1  CDR  CBDCOM 

ATTN:  TECHNICAL  LIBRARY 
BUXJ  E3330 

1  DIR  CBIAC 

BLDG  E3330,  RM  150 


1  DEPUTY  CHIEF  OF  STAFF 
OPERATIONS  AND  PLANS 
ATTN  DAMO  SW  RM  3C630 
400  ARMY  PENTAGON 
WASHINGTON  DC  20310-0400 
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NO.  OF 

COPffiS  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


1  HQDA 

DUSA  OR 

ROOM2E660 

102  ARMY  PENTAGON 

WASHINGTON  DC  20310-0102 

1  DEP  CHIEF  OF  STAFF 

OPERATIONS  AND  PLANS 
ATTN  DAMO  SW  ROOM  3C630 
400  ARMY  PENTAGON 
WASHINGTON  DC  20310-0400 

1  COMMANDER 

US  ARMY  MISSILE  COMMAND 
ATTN  AMSTA  CG 

REDSTONE  ARSENAL  AL  35898-5000 

1  DIRECTOR 

US  ARMY  BALUSTIC  MIS  DEFNS 
SYS  CMND 

ATTN  ADVNCD  TECHLGY  CTR 
PO  BOX  1500 

HUNTSVEXE  AL  35807-3801 

2  COMMANDER 

USA  STRTGC  DEFNS  CMND 
ATTN  CSSD  SL  C 
CSSD  SL  S 

HUNTSVILLE  AL  35807-3801 

1  COMMANDER 

CECOM  R&D  TEOI  UB 
ATTN  ASQNC  ELC  IS  L  R 
FORT  MONMOUTH  NJ  07703-5000 

1  COMMANDER 

US  ARMY  RESEARCH  OFFICE 
ATTN  TECH  LIB 
PO  BOX  12211 

RESEARCH  TRIANGLE  PARK  NC  27709-2211 

1  COMMANDER 

USA  LOGISTICS  MGMT  CTR 
ATTN  DEFNS  LGSTCS  STUDIES 
FORT  LEE  VA  23801 

1  COMMANDER 

US  ARMY  NGIC 
ATTN  AMXST  MC  3 
220  SEVENTH  ST  NE 
CHARLOTTESVILLE  VA  22901-53% 


1  DIRECTOR 
USARL 

ATTN  AMSRL  OP  TL 
WATERTOWN  MA  02172-0001 

DIRECTOR 

USARL 

ATTN  AMSRL  SL  E  G  MARES 

AMSRL  SL  EA  R  SHELBURNE 

AMSRL  SL  EG  J  PALOMO 

AMSRL  SL  EM  R  FLORES 

AMSRL  SL  EP  D  ALVAREZ 

AMSRL  SL  ES  T  ATHERTON 

AMSRL  SL  EV  DR  K  MORRISON 

AMSRL  SL  CA  R  SUTHERLAND 

WHITE  SANDS  MISSILE  RANGE  NM  88002-5513 

1  DIRECTOR 

USARL 

ATTN  AMSRL  EP 
M  V  GELNOVATCH 
FT  MONMOUTH  NJ  07703-5601 

1  DIRECTOR 

USARL 

ATTN  AMSRL  EP  EF 
DR  M.  DUTTA 

FT  MONMOUTH  NJ  07703-5601 

1  DIRECTOR 

USARL 

ATTN  AMSRL  SS  V  DEMONTE 
ADELPHI  MD  20783-1197 

1  DIRECTOR 

USARL 

ATTN  AMSRL  SS  IC 
DRPEMMERMAN 
ADELPHI  MD  20783-1197 

1  DIRECTOR 

USARL 

ATTN  AMSRL  SS  IC 
R  WINKLER 

ADELPHI  MD  20783-1197 

1  DIRECTOR 

USARL 

ATTN  AMSRL  WT  NH 
DR  H  E  BRANDT 
ADELPHI  MD  20783-1197 
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NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


1  COMMANDANT 

USA  CMND  AND  GENL  STAFF 
FORT  LEAVENWORTH  KS  66027 

1  DIRECTOR 

INSTITUTE  OF  DEFNS  ANLYS 
ATTN  LIBRARY 
1801  BEAUREGARD  ST 
ALEXANDRIA  VA  22311 


ARERDFEN  PROVING  GROUND 

1  DIR,  USAERDEC 
ATTN  SCBRD  RT 

1  CDR,  USACBDA 

ATTN  AMSCB  CH 

1  CDR,  USAATC 

ATTN  STECS 

77  DIR,  USARL 

ATTN  AMSRL  SL  J  SMITH 
AMSRL  SL  BA 
J  MORRISSEY 
J  HANES 
RKUNKEL 
L  ROACH 

AMSRL  SL  BG  A  YOUNG 
AMSRL  SL  BS 
DBELY 
T  KLOPCIC 
AMSRL  SL  BV 
R  SANDMEYER 
W  MERMAGEN  JR 
AMSRL  SL  C 
WJ  HUGHES 
J  SEIGH 
MAJ  GILMAN 
AMSRL  SL  CM 
DFARENWALD 
B  RUTH 
L  DAVIS 
L  DELICIO 
R  JOLLIFFE 
D  SLOOP 
RTYTUS 

R  ZUM  BRUNNEN 
E  FIORAVANTE 
MMAR 

DR  §OLN  (30  CPS) 


AMSRL  SL  CO 
D  BAYLOR 
R  PROCHAZKA 
AMSRL  SL  CS 
J  BEILFUSS 
TFLORY 
B  SMITH 
M  SMITH 
JFRANZ 
TMAK 
JNEALON 
R  PARSONS 
M  KAUFMAN 
DMANYAK 
R  POLIMADEI 
M  BUMBAUGH 
J  CAPOBIANCO 
D  DAVIS 
AMSRL  SL  I 
R  REITZ 
D  BASSETT 
A  YOUNG 
M  VOGEL 
E  PANUSKA 
D  HASKELL 
J FEENEY 
R  ZIGLER 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your  comments/answers 
to  the  items/questions  below  will  aid  us  in  our  efforts. 

1.  ARL  Report  Ntimber  ARL-TR-979 _ Date  of  Report  March  1996 _ 

2.  Date  Report  Received _ _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  puipose,  related  project,  or  other  area  of  interest  for  which  the  report 

wiU  be  used.)  _ _ _ _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of  ideas,  etc.) 


5.  Has  the  infoimation  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved,  operating  costs 
avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to 
oiganization,  technical  content,  format,  etc.) _ 


Organization 

CURRENT  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 

City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the 
Old  or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 
(DO  NOT  STAPLE) 


DEPARTMENT  OF  THE  ARMY 


OFRCtAL  BUSINESS 


BUSINESS  REPLY  MAIL 

FIRST  CLASS  PERMIT  NO  0001  ,APG,MD 


NO  POSTAGE 
NECESSARY 
IF  MAILED 
IN  THE 

UNITED  STATES 


POSTAGE  WILL  BE  PAID  BY  ADDRESSEE 


DIRECTOR 

U.S.  ARMY  RESEARCH  LABORATORY 
ATTN:  AMSRL-SL-CM 

ABERDEEN  PROVING  GROUND,  MD  21010-5423 


